using TaylorSeries
using MatrixEquations
using ForwardDiff
using Distances
using Plots
using Statistics
using NaNStatistics
using LinearAlgebra
using RootsAndPoles
using SpecialFunctions
using JLD2
using CSV
using DataFrames
using Optim
using Roots
using Random
using OrdinaryDiffEq
using QuadGK
using Distributions
using StatsBase
using DelimitedFiles
using Printf
using LaTeXStrings

using Pkg
Pkg.add(PackageSpec(url="https://github.com/StuartBenjamin/ScrewPinchAsymptoticMatching.jl"))
using ScrewPinchAsymptoticMatching
ScrewPinchAsymptoticMatching.mod_path=chop(pathof(ScrewPinchAsymptoticMatching);tail=31)

include("/Users/sbenjamin/Desktop/PHD/Cylindrical Delta Prime Widget/Screw pinch asymptotic matching in Julia/ScrewPinchAsymptoticMatching.jl/src/DataBaseGen.jl")
##########################################################################################################################################################
#Mass generating of m=2,3 for n=1: Let's get rockin
##########################################################################################################################################################

mu0 = ScrewPinchAsymptoticMatching.mu0

if true #These are a mix of ARC and SPARC parameters 
    q0 = 1.1 #Moves qstart, q still increases by the same proportion regardless of where you place it
    rs0 = 0.7 #Controls qedge to be around 3.3 (SPARC is 3.4)
    R0 = 3.0 #Major radius
    ν = 1.0  #Current peaking is fine, it'll change anyway 
                #Increase ν to increase total current by widening peak. Leaves peak magnitude unchanged
    xb = 1/rs0 #Sets minor radius to 1.0
    Bp0 = 2.5 #Sets Bt to around 12T (~ARC/SPARC) 
    β = 1e-2 #This is a SPARC number #14; #Gives a true 0-pressure Chandra Equil

    m = 2;
    ms = [2,3,4,5]
    n = 1;
    k = ScrewPinchAsymptoticMatching.k_(n,R0)
    #c0 = 1;
end

Chandra_Equil_qvar(plot_equil) = q0 -> Chandra_Equil(β,rs0,R0,ν,xb; Bp0=Bp0, q0=q0, plot_equil=plot_equil, print_mathematica_inputs=false, plotrvec = range(0.000001,xb*rs0,200))
#Furth_Equil_qvar(plot_equil) = q0 -> Furth_Equil(q0,rs0,R0,ν,xb; Bp0=Bp0, plot_equil=plot_equil, print_mathematica_inputs=false, plotrvec = range(0.000001,xb*rs0,200))
#Furth_Equil_qvar(true)(1.1)

Bp,Bt,q,dpdr,p,Jt,Jp,rb,outerp6 = Chandra_Equil_qvar(true)(1.1)

rs = find_rs(q,m,n,rb,q0,ν,rs0)
ref_equil=Equilibrium(Bp,Bt,q,dpdr,p,r->Jt(r),Jp,rb,R0,rs,rs0)
plot_equil(ref_equil)

raw_delta_prime(Bp, Bt, dpdr, k, m, r0, rs, rb, 1e-8)
Δl_Δr_calculator(Bp, Bt, dpdr, k, m, r0, rs, rb, rs0, 1, 1e-8;plot_soln=true,plot_soln_1=true)

total_plasma_current(Jt,1.0)

########################### Current generation inputs
Jtotmax = 1.1*total_plasma_current(Jt,rb)
Jtotmin = 0.9*total_plasma_current(Jt,rb)
Jtotrange=[Jtotmax, Jtotmin]
J0bounds=[1.2*Jt(0.0),0.5*Jt(0.0)]
Jedgebounds=[0.3*Jt(0.0),0.0*Jt(0.0)]
Jt_upper_bound=maximum(J0bounds)
knotmax=10
knotmin=5
dirichlet_alpha=10
maxgrad_width = rb/2
max_ref=Jt(0.0)
wobbles=3


r0 = 1e-2;
del=1e-5;
nmax=1;

"""
f_gen_equilibria(batch_size, num_clean) = gen_n_clean_equilibria(Jt_upper_bound,p,rb,batch_size,num_clean; dpdr=dpdr,Bt0=Bt(0.0), R0=R0, ideal_verbose=false, qtest=m/n, 
Jtot_range = Jtotrange, use_fine=true, use_coarse=true, monotonic=false, wobbles=wobbles, max_ref=max_ref,maxgrad_width = maxgrad_width, maxbatches=100, 
J0bounds=J0bounds, Jedgebounds=Jedgebounds, dirichlet_alpha=dirichlet_alpha, knotmax=knotmax, knotmin=knotmin);

f_gen_equilibria_6_knots(batch_size, num_clean) = gen_n_clean_equilibria(Jt_upper_bound,p,rb,batch_size,num_clean; dpdr=dpdr,Bt0=Bt(0.0), R0=R0, ideal_verbose=false, qtest=m/n, Jtot_range = Jtotrange, use_fine=true, use_coarse=true, monotonic=false, wobbles=wobbles, max_ref=max_ref,maxgrad_width = maxgrad_width, maxbatches=100, 
J0bounds=J0bounds, Jedgebounds=Jedgebounds, dirichlet_alpha=dirichlet_alpha, knotmax=knotmax, knotmin=6);

f_gen_equilibria_7_knots(batch_size, num_clean) = gen_n_clean_equilibria(Jt_upper_bound,p,rb,batch_size,num_clean; dpdr=dpdr,Bt0=Bt(0.0), R0=R0, ideal_verbose=false, qtest=m/n, Jtot_range = Jtotrange, use_fine=true, use_coarse=true, monotonic=false, wobbles=wobbles, max_ref=max_ref,maxgrad_width = maxgrad_width, maxbatches=100, 
J0bounds=J0bounds, Jedgebounds=Jedgebounds, dirichlet_alpha=dirichlet_alpha, knotmax=knotmax, knotmin=7);

f_gen_equilibria_8_knots(batch_size, num_clean) = gen_n_clean_equilibria(Jt_upper_bound,p,rb,batch_size,num_clean; dpdr=dpdr,Bt0=Bt(0.0), R0=R0, ideal_verbose=false, qtest=m/n, Jtot_range = Jtotrange, use_fine=true, use_coarse=true, monotonic=false, wobbles=wobbles, max_ref=max_ref,maxgrad_width = maxgrad_width, maxbatches=100, 
J0bounds=J0bounds, Jedgebounds=Jedgebounds, dirichlet_alpha=dirichlet_alpha, knotmax=knotmax, knotmin=8);
"""

cd("/Users/sbenjamin/Desktop/PHD/ScrewPinchDataBase")
path="/Users/sbenjamin/Desktop/PHD/ScrewPinchDataBase"

batch_size=100000
batch_size6=4*batch_size
batch_size7=100*batch_size
batch_size8=1000*batch_size
num_clean=1000
loop_big_batch=5

#gen_equil_run_Δl_Δr(f_gen_equilibria_6_knots, batch_size6, num_clean, ms, n, r0, nmax, del; filename = "n2_study_fin_pressure_6knot_1.jld2", loop_big_batch=loop_big_batch, path=path)
#gen_equil_run_Δl_Δr(f_gen_equilibria_6_knots, batch_size6, num_clean, ms, n, r0, nmax, del; filename = "n2_study_fin_pressure_6knot_2.jld2", loop_big_batch=loop_big_batch, path=path)
#gen_equil_run_Δl_Δr(f_gen_equilibria_7_knots, batch_size7, num_clean, ms, n, r0, nmax, del; filename = "n2_study_fin_pressure_7knot_1.jld2", loop_big_batch=loop_big_batch, path=path)
#gen_equil_run_Δl_Δr(f_gen_equilibria_7_knots, batch_size7, num_clean, ms, n, r0, nmax, del; filename = "n2_study_fin_pressure_7knot_2.jld2", loop_big_batch=loop_big_batch, path=path)
#gen_equil_run_Δl_Δr(f_gen_equilibria_8_knots, batch_size8, num_clean, ms, n, r0, nmax, del; filename = "n2_study_fin_pressure_8knot_1.jld2", loop_big_batch=loop_big_batch, path=path)
#gen_equil_run_Δl_Δr(f_gen_equilibria_6_knots, batch_size6, num_clean, ms, n, r0, nmax, del; filename = "n2_study_fin_pressure_6knot_3.jld2", loop_big_batch=loop_big_batch, path=path)
#gen_equil_run_Δl_Δr(f_gen_equilibria_7_knots, batch_size7, num_clean, ms, n, r0, nmax, del; filename = "n2_study_fin_pressure_7knot_3.jld2", loop_big_batch=loop_big_batch, path=path)
#gen_equil_run_Δl_Δr(f_gen_equilibria_8_knots, batch_size8, num_clean, ms, n, r0, nmax, del; filename = "n2_study_fin_pressure_8knot_2.jld2", loop_big_batch=loop_big_batch, path=path)
#gen_equil_run_Δl_Δr(f_gen_equilibria_8_knots, batch_size8, num_clean, ms, n, r0, nmax, del; filename = "n2_study_fin_pressure_8knot_3.jld2", loop_big_batch=loop_big_batch, path=path)

#gen_equil_run_Δl_Δr(f_gen_equilibria, batch_size, num_clean, ms, n, r0, nmax, del; filename = "n2_study_fin_pressure1.jld2", loop_big_batch=loop_big_batch, path=path)
#gen_equil_run_Δl_Δr(f_gen_equilibria, batch_size, num_clean, ms, n, r0, nmax, del; filename = "n2_study_fin_pressure2.jld2", loop_big_batch=loop_big_batch, path=path)
#gen_equil_run_Δl_Δr(f_gen_equilibria, batch_size, num_clean, ms, n, r0, nmax, del; filename = "n2_study_fin_pressure3.jld2", loop_big_batch=loop_big_batch, path=path)
#gen_equil_run_Δl_Δr(f_gen_equilibria, batch_size, num_clean, ms, n, r0, nmax, del; filename = "n2_study_fin_pressure4.jld2", loop_big_batch=loop_big_batch, path=path)
#gen_equil_run_Δl_Δr(f_gen_equilibria, batch_size, num_clean, ms, n, r0, nmax, del; filename = "n2_study_fin_pressure5.jld2", loop_big_batch=loop_big_batch, path=path)
#gen_equil_run_Δl_Δr(f_gen_equilibria, batch_size, num_clean, ms, n, r0, nmax, del; filename = "n2_study_fin_pressure6.jld2", loop_big_batch=loop_big_batch, path=path)
#gen_equil_run_Δl_Δr(f_gen_equilibria, batch_size, num_clean, ms, n, r0, nmax, del; filename = "n2_study_fin_pressure7.jld2", loop_big_batch=loop_big_batch, path=path)
#gen_equil_run_Δl_Δr(f_gen_equilibria, batch_size, num_clean, ms, n, r0, nmax, del; filename = "n2_study_fin_pressure8.jld2", loop_big_batch=loop_big_batch, path=path)
#gen_equil_run_Δl_Δr(f_gen_equilibria, batch_size, num_clean, ms, n, r0, nmax, del; filename = "n2_study_fin_pressure9.jld2", loop_big_batch=loop_big_batch, path=path)
#gen_equil_run_Δl_Δr(f_gen_equilibria, batch_size, num_clean, ms, n, r0, nmax, del; filename = "n2_study_fin_pressure10.jld2", loop_big_batch=loop_big_batch, path=path)

########################### Analysis
cd("/Users/sbenjamin/Desktop/PHD/ScrewPinchDataBase")
resistive_equils_n2 = ResistiveEquilibrium[]
@load "n2_study_fin_pressure1.jld2" resistive_equils_store equils_and_inds
for i in resistive_equils_store
    append!(resistive_equils_n2,i)
end
@load "n2_study_fin_pressure2.jld2" resistive_equils_store equils_and_inds
for i in resistive_equils_store
    append!(resistive_equils_n2,i)
end
@load "n2_study_fin_pressure3.jld2" resistive_equils_store equils_and_inds
for i in resistive_equils_store
    append!(resistive_equils_n2,i)
end
@load "n2_study_fin_pressure4.jld2" resistive_equils_store equils_and_inds
for i in resistive_equils_store
    append!(resistive_equils_n2,i)
end
@load "n2_study_fin_pressure5.jld2" resistive_equils_store equils_and_inds
for i in resistive_equils_store
    append!(resistive_equils_n2,i)
end
@load "n2_study_fin_pressure6.jld2" resistive_equils_store equils_and_inds
for i in resistive_equils_store
    append!(resistive_equils_n2,i)
end
@load "n2_study_fin_pressure7.jld2" resistive_equils_store equils_and_inds
for i in resistive_equils_store
    append!(resistive_equils_n2,i)
end
@load "n2_study_fin_pressure_6knot_1.jld2" resistive_equils_store equils_and_inds
for i in resistive_equils_store
    append!(resistive_equils_n2,i)
end
@load "n2_study_fin_pressure_6knot_2.jld2" resistive_equils_store equils_and_inds
for i in resistive_equils_store
    append!(resistive_equils_n2,i)
end
@load "n2_study_fin_pressure_7knot_1.jld2" resistive_equils_store equils_and_inds
for i in resistive_equils_store
    append!(resistive_equils_n2,i)
end
@load "n2_study_fin_pressure_7knot_2.jld2" resistive_equils_store equils_and_inds
for i in resistive_equils_store
    append!(resistive_equils_n2,i)
end

length(resistive_equils_n2)

########################### Re-doing some analyses with higher m, n.
###########################

"""

##resistive_equils43_filtered = filter(i->q_one_rs(i.equilibrium, 4/3),resistive_equils_n2)
##@save "43_mid_step_filtered_resistive_equils.jld2" resistive_equils43_filtered
@load "43_mid_step_filtered_resistive_equils.jld2" resistive_equils43_filtered

equils_43 = [i.equilibrium for i in resistive_equils43_filtered]
Δprimes_43 = [i.Δprimes for i in resistive_equils43_filtered]


########################### 32 check
equils_43_32rss = [find_rs(i.q,3,2,i.Jt.xs[end];verbose=false) for i in collect(equils_43[1:10000])]
sorted_inds_32_rs = sortperm(equils_43_32rss)

equils_43_32_10000 = collect(equils_43[1:10000])

equils_43_32_10000[sorted_inds_32_rs[1:10]]

resistive_equils32_rssorted, inds32_rssorted = run_Δl_Δr_calculator(equils_43_32_10000[sorted_inds_32_rs[1:100]], [3], 2, r0, nmax, del; report_err=true)


deltaprimes_mn32_rssorted = [i.Δprimes[1].Δprime for i in resistive_equils32_rssorted]
rss_mn32_rssorted = [i.Δprimes[1].rs for i in resistive_equils32_rssorted]
scatter(rss_mn32_rssorted,deltaprimes_mn32_rssorted)

########################### 43 check
equils_43rss = [find_rs(i.q,4,3,i.Jt.xs[end];verbose=false) for i in collect(equils_43[1:10000])]
sorted_inds_43rs = sortperm(equils_43_32rss)

equils_43_10000 = collect(equils_43[1:10000])

equils_43_10000[sorted_inds_43rs[1:10]]

resistive_equils43_rssorted, inds43_rssorted = run_Δl_Δr_calculator(equils_43_10000[sorted_inds_43rs[1:10]], [4], 3, r0, nmax, del; report_err=true)

deltaprimes_mn43_rssorted = [i.Δprimes[1].Δprime for i in resistive_equils43_rssorted]
rss_mn43_rssorted = [i.Δprimes[1].rs for i in resistive_equils43_rssorted]
scatter(rss_mn43_rssorted,deltaprimes_mn43_rssorted)

########################### Junk:

for i in Δprimes_43[1:10000]
    rs = find_rs(i.q,3,2,i.Jt.xs[end];verbose=false)
    push!()

n = 2;
k = ScrewPinchAsymptoticMatching.k_(n,R0)
resistive_equils32_1000, inds32_1000 = run_Δl_Δr_calculator(equils_43[1:1000], [3], 2, r0, nmax, del; report_err=true)

deltaprimes_mn32 = [i.Δprimes[1].Δprime for i in resistive_equils32_1000]
rss_mn32 = [i.Δprimes[1].rs for i in resistive_equils32_1000]
sorted_inds_mn3 = sortperm(deltaprimes_mn32)

histogram2d(rss_mn32,deltaprimes_mn32; xlabel="3/2 rational surface location (m)",ylabel = "Δ' 3/2")
corspearman(rss_mn32,deltaprimes_mn32)
cor(rss_mn32,deltaprimes_mn32)

scatter(rss_mn32,deltaprimes_mn32)

run_Δl_Δr_calculator(equils_43[1:10], 2, 1, r0, nmax, del; report_err=false)

equils_n2 = [i.equilibrium for i in resistive_equils_n2]
current_Δprimes = [i.Δprimes for i in resistive_equils_n2]

resistive_equils54, inds54 = run_Δl_Δr_calculator(equils_n2, 5, 4, r0, nmax, del; report_err=false)

########################### Analyses with higher m, n's ENDED
###########################
"""

########################### Removing scenarios without n=3 rs (just the tail end of Jtot distribution -> largest currents cause qedge to drop down)
###########################
knotnums = [length(i.equilibrium.Jt.xs) for i in resistive_equils_n2]
deltaprimes_n2 = [i.Δprimes[1].Δprime for i in resistive_equils_n2]
#deltaprimes_n2 = filter(i -> resistive_equils_n2[i].Δprimes[1].m==2, 1:length(deltaprimes_n2))
deltaprimes_n3_inds = filter(i -> (length(resistive_equils_n2[i].Δprimes)>1 && resistive_equils_n2[i].Δprimes[2].m==3), 1:length(resistive_equils_n2))
deltaprimes_no_n3_inds = filter(i -> (length(resistive_equils_n2[i].Δprimes)==1), 1:length(resistive_equils_n2))
deltaprimes_n3 = [i.Δprimes[2].Δprime for i in resistive_equils_n2[deltaprimes_n3_inds]]

Jtots = [total_plasma_current(i.equilibrium.Jt,i.equilibrium.rb) for i in resistive_equils_n2]
mean(Jtots[deltaprimes_n3_inds])
mean(Jtots[deltaprimes_no_n3_inds])
plot_current_profiles(resistive_equils_n2[deltaprimes_no_n3_inds])
plot_current_profiles(resistive_equils_n2[deltaprimes_n3_inds][1:300])

resistive_equils_n2_n3 = filter(i -> (length(i.Δprimes)>1), resistive_equils_n2)
deltaprimes_n2 = [i.Δprimes[1].Δprime for i in resistive_equils_n2_n3]
deltaprimes_n3 = [i.Δprimes[2].Δprime for i in resistive_equils_n2_n3]
deltaprime_diff = deltaprimes_n2-deltaprimes_n3

########################### Removing n2 delta' outliers (they're at very small rs, integrator often broken at these rs)

rss_n2 = [i.Δprimes[1].rs for i in resistive_equils_n2_n3]
sorted_inds_n2 = sortperm(deltaprimes_n2)
sorted_inds_n3 = sortperm(deltaprimes_n3)

histogram2d(rss_n2,deltaprimes_n2; xlabel="2/1 rational surface location (m)",ylabel = "Δ' 2/1")
corspearman(rss_n2,deltaprimes_n2)
cor(rss_n2,deltaprimes_n2)


deltaprimes_n2[sorted_inds_n2[1:90]]
maximum(rss_n2[sorted_inds_n2[1:90]])

deltaprimes_n2[sorted_inds_n2[80:90]]
rss_n2[sorted_inds_n2[80:90]]


iend = 87
deltaprimes_n2[sorted_inds_n2[iend]]
maximum(rss_n2[sorted_inds_n2[1:iend]])

scatter(rss_n2[sorted_inds_n2[1:iend]],deltaprimes_n2[sorted_inds_n2[1:iend]]; xlabel="2/1 rational surface location (m)",ylabel = "Δ' 2/1")
corspearman(rss_n2[sorted_inds_n2[1:iend]],deltaprimes_n2[sorted_inds_n2[1:iend]])
cor(rss_n2[sorted_inds_n2[1:iend]],deltaprimes_n2[sorted_inds_n2[1:iend]])

resistive_equils_n2_n3_clean = filter(i -> (i.Δprimes[1].Δprime>deltaprimes_n2[sorted_inds_n2[iend]]), resistive_equils_n2_n3)
deltaprimes_n2 = [i.Δprimes[1].Δprime for i in resistive_equils_n2_n3_clean]
deltaprimes_n3 = [i.Δprimes[2].Δprime for i in resistive_equils_n2_n3_clean]
deltaprime_diff = deltaprimes_n2-deltaprimes_n3

########################### Cases where n2 is stable and n3 is unstable (rare)
negm2posm3_inds = filter(i -> deltaprimes_n2[i]<0.0&& deltaprimes_n3[i]>0.0, 1:length(resistive_equils_n2_n3_clean))
negm2posm3_equils = filter(i -> i.Δprimes[1].Δprime<0.0 && i.Δprimes[2].Δprime>0.0, resistive_equils_n2_n3_clean)
plot_current_profiles(negm2posm3_equils)

########################### general analysis, n2, n3 histograms
histogram(deltaprimes_n2)
histogram(filter(x->x>-200,deltaprimes_n3))
n_3_inds = filter(i -> (deltaprimes_n3[i]>-200), 1:length(deltaprimes_n3))

Jts = [i.equilibrium.Jt for i in resistive_equils_n2_n3_clean]
Jps = [i.equilibrium.Jp for i in resistive_equils_n2_n3_clean]
rss_n2 = [i.Δprimes[1].rs for i in resistive_equils_n2_n3_clean]
rss_n3 = [i.Δprimes[2].rs for i in resistive_equils_n2_n3_clean]
Jtots = [total_plasma_current(i.equilibrium.Jt,i.equilibrium.rb) for i in resistive_equils_n2_n3_clean]
Jtgrad_at_rs_n2 = [ScrewPinchAsymptoticMatching.gradient(Jts[i],rss_n2[i]) for i in 1:length(Jts)]
Jtgrad_at_rs_n3 = [ScrewPinchAsymptoticMatching.gradient(Jts[i],rss_n3[i]) for i in 1:length(Jts)]
Jpgrad_at_rs_n2 = [ForwardDiff.derivative(Jps[i],rss_n2[i]) for i in 1:length(Jps)]
Jpgrad_at_rs_n3 = [ForwardDiff.derivative(Jps[i],rss_n3[i]) for i in 1:length(Jps)]
Jt_at_rs_n2 = [Jts[i](rss_n2[i]) for i in 1:length(Jts)]
Jt_at_rs_n3 = [Jts[i](rss_n3[i]) for i in 1:length(Jts)]
shear_at_rs_n2 = [ForwardDiff.derivative(resistive_equils_n2_n3_clean[i].equilibrium.q,rss_n2[i]) for i in 1:length(Jts)]
shear_at_rs_n3 = [ForwardDiff.derivative(resistive_equils_n2_n3_clean[i].equilibrium.q,rss_n3[i]) for i in 1:length(Jts)]
p_at_rs_n2 = [mu0*resistive_equils_n2_n3_clean[i].equilibrium.p(rss_n2[i]) for i in 1:length(resistive_equils_n2_n3_clean)]
p_at_rs_n3 = [mu0*resistive_equils_n2_n3_clean[i].equilibrium.p(rss_n3[i]) for i in 1:length(resistive_equils_n2_n3_clean)]
dpdr_at_rs_n2 = [resistive_equils_n2_n3_clean[i].equilibrium.dpdr(rss_n2[i]) for i in 1:length(resistive_equils_n2_n3_clean)]
dpdr_at_rs_n3 = [resistive_equils_n2_n3_clean[i].equilibrium.dpdr(rss_n3[i]) for i in 1:length(resistive_equils_n2_n3_clean)]
Drs_n2 = [Dr(i)[1] for i in resistive_equils_n2_n3_clean]
Drs_n3 = [Dr(i)[2] for i in resistive_equils_n2_n3_clean]

xs_s_raw = [i.equilibrium.Jt.xs for i in resistive_equils_n2_n3_clean]
knotnum_lengths = length.(xs_s_raw)

rando = rand(1:length(resistive_equils_n2_n3_clean),340)
plot_current_profiles(resistive_equils_n2_n3_clean[rando])
plot_q_profiles(resistive_equils_n2_n3_clean[rando])

histogram(Drs_n2; label=false, title=false, thickness_scaling=1.5)
histogram(knotnum_lengths; label=false, title="Distribution of radial knot number in database", thickness_scaling=1., color=:black, thickness_scaling=1.5)
vals_hist = 5:9
counts = zeros(Int,5)
for i in knotnum_lengths
    for (io,j) in enumerate(vals_hist)
        if i==j
            counts[io]+=1
        end
    end
end
counts_norm = counts./sum(counts)
bar(vals_hist,counts; label=false, title="Distribution of radial knot number in database", color=:black, thickness_scaling=1.5)#,margin=11Plots.mm)
plot_current_profiles(resistive_equils_n2_n3_clean[78000];plot_knots=true)
plot_equil(resistive_equils_n2_n3_clean[78000])
plot_equil_short(resistive_equils_n2_n3_clean[78000].equilibrium)

histogram(xs_s; label=false, title="Radial distribution of knots in database", thickness_scaling=1.3, bins = 200, color=:black)
save()
xs_7s_raw = filter(i -> length(i)>6, xs_s_raw)
xs_s = Number[]
xs_7 = Number[]
for i in xs_s_raw
    append!(xs_s, i[2:end-1])
end
for i in xs_7s_raw
    append!(xs_7, i[2:end-1])
end
histogram(xs_s)
histogram(xs_7)

histogram(rss_n2; xlabel="radius (m)", label=false, xlims=(0.0,1.0), bins=200, title="Radial distribution of 2/1 rational surface", color=:black, thickness_scaling=1.4)
histogram(rss_n3; xlabel="radius (m)", label=false, xlims=(0.0,1.0), bins=200, title="Radial distribution of 3/1 rational surface", color=:black, thickness_scaling=1.4)

histogram([rss_n2,rss_n3]; xlabel="radius (m)", label=false, xlims=(0.0,1.0), bins=100, title="Radial distribution of 2/1 rational surface",  thickness_scaling=1.2)

histogram([rss_n3,rss_n2]; xlabel="radius (m)", label=false, xlims=(0.0,1.0), bins=100, title="Radial distribution of 2/1 rational surface",  thickness_scaling=1.2)


!histogram(rss_n3; xlabel="radius (m)", label=false, xlims=(0.0,1.0), bins=200, title="Radial distribution of 3/1 rational surface",  thickness_scaling=1.2)


histogram2d(rss_n2[negm2posm3_inds],rss_n3[negm2posm3_inds],xlabel="2/1 rational surface location (m)",ylabel = "3/1 rational surface location (m)")

histogram(Jtots./1e6; xlabel="\$I\$ (MA)", label=false, title="Distribution of total toroidal current", bins=200, color=:black, thickness_scaling=1.4)
histogram(knotnums)

########################### Important correlation: rs approaching edge of plasma

histogram2d(rss_n2,deltaprimes_n2; xlabel="2/1 rational surface location (m)",ylabel = "Δ' 2/1")
cor(rss_n2,deltaprimes_n2)
corspearman(rss_n2,deltaprimes_n2)
histogram2d(rss_n3[n_3_inds],deltaprimes_n3[n_3_inds]; xlabel="3/1 rational surface location (m)",ylabel = "Δ' 3/1") 
cor(rss_n3[n_3_inds],deltaprimes_n3[n_3_inds])
corspearman(rss_n3[n_3_inds],deltaprimes_n3[n_3_inds])
#^DFINITELY A CONSEQUENCE OF ENFORCING WE HAVE AN INTERNAL MODE (perturbation = 0 at edge)

########################### n2, n3 extreme Cases
sorted_inds_n2 = sortperm(deltaprimes_n2)
sorted_inds_n3 = sortperm(deltaprimes_n3)

deltaprimes_n2[sorted_inds_n2[1:10]]
deltaprimes_n3[sorted_inds_n3[1:10]]

deltaprimes_n2[sorted_inds_n2[end-10:end]]
deltaprimes_n3[sorted_inds_n3[end-10:end]]

histogram(deltaprimes_n2)

plot_current_profiles(resistive_equils_n2_n3_clean[sorted_inds_n2[end-10:end]];m_ind=1) #THESE ARE EXTREMELY FLAT CURRENT PROFILES pushing rs_n2 towards 0 -> not interesting
plot_equil(resistive_equils_n2_n3_clean[sorted_inds_n2[end-10:end]];case=1, Jt_ref = 1.1*Jt_upper_bound)
plot_current_profiles(resistive_equils_n2_n3_clean[sorted_inds_n2[1:10]];m_ind=1) #THESE ARE the rs_n2 on the hill -> INTERESTING
plot_current_profiles(resistive_equils_n2_n3_clean[sorted_inds_n3[1:10]];m_ind=2) #THESE rs_n3 towards rb stabilising mode -> not interesting
plot_current_profiles(resistive_equils_n2_n3_clean[sorted_inds_n3[end-10:end]];m_ind=2)  #THESE ARE the rs_n3 after long flat period, unstable n3 -> INTERESTING



########################### Cases where n_2 rs values aren't in super low-rs region:
histogram2d(rss_n2,deltaprimes_n2; xlabel="2/1 rational surface location (m)",ylabel = "Δ' 2/1")
cor(rss_n2,deltaprimes_n2)
histogram2d(rss_n3[n_3_inds],deltaprimes_n3[n_3_inds]; xlabel="3/1 rational surface location (m)",ylabel = "Δ' 3/1")
cor(rss_n3[n_3_inds],deltaprimes_n3[n_3_inds])

non_rs_inds = filter(i -> resistive_equils_n2_n3[i].Δprimes[1].rs>0.3, 1:length(resistive_equils_n2_n3_clean)) #getting rid of any correlation with position for n2
non_rs_equil = filter(i -> i.Δprimes[1].rs>0.3, resistive_equils_n2_n3_clean)

########################### First big thing: local gradient correlates with Delta' for 2/1 (evidence it might be absolute total local current gradient)
deltaprimes_n2_nrs= [i.Δprimes[1].Δprime for i in non_rs_equil]
deltaprimes_n3_nrs = [i.Δprimes[2].Δprime for i in non_rs_equil]
sorted_inds_n2_nrs = sortperm(deltaprimes_n2_nrs)
sorted_inds_n3_nrs = sortperm(deltaprimes_n3_nrs)

Jts_nrs = [i.equilibrium.Jt for i in non_rs_equil]
Jps_nrs = [i.equilibrium.Jp for i in non_rs_equil]
rss_n2_nrs = [i.Δprimes[1].rs for i in non_rs_equil]
rss_n3_nrs = [i.Δprimes[2].rs for i in non_rs_equil]
Jtots_nrs = [total_plasma_current(i.equilibrium.Jt,i.equilibrium.rb) for i in non_rs_equil]

histogram(deltaprimes_n2_nrs;xlabel="Δ' 2/1", label=false, bins=1000, title = "Δ' 2/1 distribution")

histogram2d(rss_n2_nrs,deltaprimes_n2_nrs; xlabel="2/1 rational surface location (m)",ylabel = "Δ' 2/1")
cor(rss_n2_nrs,deltaprimes_n2_nrs)
corspearman(rss_n2_nrs,deltaprimes_n2_nrs)
n_3_inds_temp = filter(i -> (deltaprimes_n3_nrs[i]>-200), 1:length(deltaprimes_n3_nrs))
histogram2d(rss_n3_nrs[n_3_inds_temp],deltaprimes_n3_nrs[n_3_inds_temp]; xlabel="3/1 rational surface location (m)",ylabel = "Δ' 3/1")
cor(rss_n3_nrs[n_3_inds_temp],deltaprimes_n3_nrs[n_3_inds_temp])
corspearman(rss_n3_nrs[n_3_inds_temp],deltaprimes_n3_nrs[n_3_inds_temp])

Jtgrad_at_rs_n2_nrs = [ScrewPinchAsymptoticMatching.gradient(Jts_nrs[i],rss_n2_nrs[i]) for i in 1:length(non_rs_equil)]
Jtgrad_at_rs_n3_nrs = [ScrewPinchAsymptoticMatching.gradient(Jts_nrs[i],rss_n3_nrs[i]) for i in 1:length(non_rs_equil)]
Jpgrad_at_rs_n2_nrs = [ForwardDiff.derivative(Jps_nrs[i],rss_n2_nrs[i]) for i in 1:length(non_rs_equil)]
Jpgrad_at_rs_n3_nrs = [ForwardDiff.derivative(Jps_nrs[i],rss_n3_nrs[i]) for i in 1:length(non_rs_equil)]

plot_current_profiles(non_rs_equil[sorted_inds_n2_nrs[1:5]];m_ind=1)
deltaprimes_n2_nrs[sorted_inds_n2_nrs[1:5]]
plot_current_profiles(non_rs_equil[sorted_inds_n2_nrs[1:10]];m_ind=1,color=:black)

plot_current_profiles(non_rs_equil[sorted_inds_n3_nrs[1:10]];m_ind=2)
deltaprimes_n3_nrs[sorted_inds_n3_nrs[1:10]]

########################### Terrible stable Psi: 
great_equils = collect([i.equilibrium for i in non_rs_equil[sorted_inds_n2_nrs[1:10]]])
#run_Δl_Δr_calculator(teriible_equils[1:2], [2], 1, r0, nmax, del; report_err=true, plot_soln_equil=true)

i=great_equils[3]
rsi = find_rs(i.q,2,1,i.Jt.xs[end];verbose=false)
n = 1;
k = ScrewPinchAsymptoticMatching.k_(n,R0)
Δl,Δr,del2 = Δl_Δr_calculator(i.Bp, i.Bt, i.dpdr, k, 2, r0, rsi, i.rb, i.rs0, nmax, del; integrator_reltol=1e-20, plot_solution=true, plot_soln = false)


########################### Terrible: 
plot_current_profiles(non_rs_equil[sorted_inds_n2_nrs[end-10:end]];m_ind=1) 
deltaprimes_n2_nrs[sorted_inds_n2_nrs[end-10:end]]

plot_current_profiles(non_rs_equil[sorted_inds_n2_nrs[end-10:end]];m_ind=1,color=:black)

########################### Terrible plotting Psi: 
teriible_equils = collect([i.equilibrium for i in non_rs_equil[sorted_inds_n2_nrs[end-3:end]]])
#run_Δl_Δr_calculator(teriible_equils[1:2], [2], 1, r0, nmax, del; report_err=true, plot_soln_equil=true)

i=teriible_equils[1]
rsi = find_rs(i.q,2,1,i.Jt.xs[end];verbose=false)
n = 1;
k = ScrewPinchAsymptoticMatching.k_(n,R0)
Δl,Δr,del2 = Δl_Δr_calculator(i.Bp, i.Bt, i.dpdr, k, 2, r0, rsi, i.rb, i.rs0, nmax, del; integrator_reltol=1e-20, plot_solution=true, plot_soln = false)

###GENERAL STABLE
i=8700
plot_current_profiles(non_rs_equil[sorted_inds_n2_nrs[i-100:i]];m_ind=1) 
deltaprimes_n2_nrs[sorted_inds_n2_nrs[i-100:i]]
plot_current_profiles(non_rs_equil[sorted_inds_n2_nrs[i-50:i]];m_ind=1,color=:black)

###GENERAL UNSTABLE
i=95000
plot_current_profiles(non_rs_equil[sorted_inds_n2_nrs[i-100:i]];m_ind=1) 
deltaprimes_n2_nrs[sorted_inds_n2_nrs[i-100:i]]
plot_current_profiles(non_rs_equil[sorted_inds_n2_nrs[i-100:i]];m_ind=1,color=:black)


#Jt n2 big correlation
histogram2d(Jtgrad_at_rs_n2_nrs./1e6,deltaprimes_n2_nrs; xlabel="\$J_t^'(r_s) (\\text{MA/m}^{3}) \\textrm{at } \\textrm{r_s}\$ of 2/1",ylabel = "Δ' of 2/1")
cor(Jtgrad_at_rs_n2_nrs,deltaprimes_n2_nrs)
corspearman(Jtgrad_at_rs_n2_nrs,deltaprimes_n2_nrs)

                ###########NEED TO REMOVE WELLS AND HILLS AND LOOK AGAIN!!!!!
    






#Jp n2 smaller correlation, might come from Jt (since they're all correlated)
histogram2d(Jpgrad_at_rs_n2_nrs./1e6,deltaprimes_n2_nrs; xlabel="Jp' (\$\\textrm{MA/m}^{3}\$) at rs of 2/1",ylabel = "Δ' of 2/1")
cor(Jpgrad_at_rs_n2_nrs,deltaprimes_n2_nrs)
corspearman(Jpgrad_at_rs_n2_nrs,deltaprimes_n2_nrs)

########################### Does local gradient correlate with Delta' for 3/1? Yes but other way (is it absolute total local current gradient?)
non_rs__n3_inds = filter(i -> resistive_equils_n2_n3[i].Δprimes[1].rs>0.3 && resistive_equils_n2_n3[i].Δprimes[2].rs<0.7, 1:length(resistive_equils_n2_n3_clean)) #getting rid of any correlation with position for n2
non_rs__n3_equil = filter(i -> i.Δprimes[1].rs>0.3 && i.Δprimes[2].rs<0.7, resistive_equils_n2_n3_clean)

deltaprimes_n2_nrs_n3= [i.Δprimes[1].Δprime for i in non_rs__n3_equil]
deltaprimes_n3_nrs_n3 = [i.Δprimes[2].Δprime for i in non_rs__n3_equil]
sorted_inds_n2_nrs_n3 = sortperm(deltaprimes_n2_nrs_n3)
sorted_inds_n3_nrs_n3 = sortperm(deltaprimes_n3_nrs_n3)

Jts_nrs_n3 = [i.equilibrium.Jt for i in non_rs__n3_equil]
Jps_nrs_n3 = [i.equilibrium.Jp for i in non_rs__n3_equil]
rss_n2_nrs_n3 = [i.Δprimes[1].rs for i in non_rs__n3_equil]
rss_n3_nrs_n3 = [i.Δprimes[2].rs for i in non_rs__n3_equil]
Jtots_nrs_n3 = [total_plasma_current(i.equilibrium.Jt,i.equilibrium.rb) for i in non_rs__n3_equil]

#There's a small number of equilibria where rs_n3 is small enough that the correlation with the edge goes away
mean(Jtots_nrs)
mean(Jtots_nrs_n3)
histogram(Jtots_nrs)
histogram!(Jtots_nrs_n3)
#However, there is strong trend in likelihood of rs being towards outer edge, and likelihood of total current being small
histogram(rss_n3_nrs_n3)
#these effects may be combining to confound the 'opposite trend in local gradient correlating with Delta'. However we can't be sure. Results are as they stand.

histogram2d(rss_n2_nrs_n3,deltaprimes_n2_nrs_n3; xlabel="2/1 rational surface location (m)",ylabel = "Δ' 2/1")
cor(rss_n2_nrs_n3,deltaprimes_n2_nrs_n3)
corspearman(rss_n2_nrs_n3,deltaprimes_n2_nrs_n3)
histogram2d(rss_n3_nrs_n3,deltaprimes_n3_nrs_n3; xlabel="3/1 rational surface location (m)",ylabel = "Δ' 3/1")
cor(rss_n3_nrs_n3,deltaprimes_n3_nrs_n3)
corspearman(rss_n3_nrs_n3,deltaprimes_n3_nrs_n3)

Jtgrad_at_rs_n2_nrs_n3 = [ScrewPinchAsymptoticMatching.gradient(Jts_nrs_n3[i],rss_n2_nrs_n3[i]) for i in 1:length(non_rs__n3_equil)]
Jtgrad_at_rs_n3_nrs_n3 = [ScrewPinchAsymptoticMatching.gradient(Jts_nrs_n3[i],rss_n3_nrs_n3[i]) for i in 1:length(non_rs__n3_equil)]
Jpgrad_at_rs_n2_nrs_n3 = [ForwardDiff.derivative(Jps_nrs_n3[i],rss_n2_nrs_n3[i]) for i in 1:length(non_rs__n3_equil)]
Jpgrad_at_rs_n3_nrs_n3 = [ForwardDiff.derivative(Jps_nrs_n3[i],rss_n3_nrs_n3[i]) for i in 1:length(non_rs__n3_equil)]

plot_current_profiles(non_rs__n3_equil[sorted_inds_n3_nrs_n3[1:10]];m_ind=2)

###Recheck n2s: Still there but weaker
histogram2d(Jtgrad_at_rs_n2_nrs_n3,deltaprimes_n2_nrs_n3; xlabel="∇Jt (\$\\textrm{A/m}^{2}\$) at rs of 2/1",ylabel = "Δ' of 2/1")
cor(Jtgrad_at_rs_n2_nrs_n3,deltaprimes_n2_nrs_n3)
corspearman(Jtgrad_at_rs_n2_nrs_n3,deltaprimes_n2_nrs_n3)

histogram2d(Jpgrad_at_rs_n2_nrs_n3,deltaprimes_n2_nrs_n3; xlabel="∇Jp (\$\\textrm{A/m}^{2}\$) at rs of 2/1",ylabel = "Δ' of 2/1")
cor(Jpgrad_at_rs_n2_nrs_n3,deltaprimes_n2_nrs_n3)
corspearman(Jpgrad_at_rs_n2_nrs_n3,deltaprimes_n2_nrs_n3)

###Recheck n3s: NOW STRONG
histogram2d(Jtgrad_at_rs_n3_nrs_n3,deltaprimes_n3_nrs_n3; xlabel="∇Jt (\$\\textrm{A/m}^{2}\$) at rs of 2/1",ylabel = "Δ' of 2/1")
cor(Jtgrad_at_rs_n3_nrs_n3,deltaprimes_n3_nrs_n3)
corspearman(Jtgrad_at_rs_n3_nrs_n3,deltaprimes_n3_nrs_n3)
##
histogram2d(Jpgrad_at_rs_n3_nrs_n3,deltaprimes_n3_nrs_n3; xlabel="∇Jp (\$\\textrm{A/m}^{2}\$) at rs of 2/1",ylabel = "Δ' of 2/1")
cor(Jpgrad_at_rs_n3_nrs_n3,deltaprimes_n3_nrs_n3)
corspearman(Jpgrad_at_rs_n3_nrs_n3,deltaprimes_n3_nrs_n3)

plot_current_profiles(non_rs__n3_equil[sorted_inds_n3_nrs_n3[1:10]];m_ind=2) #Good = on a hill
plot_current_profiles(non_rs__n3_equil[sorted_inds_n3_nrs_n3[1:10]];m_ind=2) #Terrible = in a stepp well, or on the side of a steep well

########################### Checking Absolute current gradient: 

#Jt n2 big correlation is unchanged
histogram2d(abs.(Jtgrad_at_rs_n2_nrs),deltaprimes_n2_nrs; xlabel="∇Jt (\$\\textrm{A/m}^{2}\$) at rs of 2/1",ylabel = "Δ' of 2/1")
cor(abs.(Jtgrad_at_rs_n2_nrs),deltaprimes_n2_nrs)
corspearman(abs.(Jtgrad_at_rs_n2_nrs),deltaprimes_n2_nrs)

#Jp n2 sligthly worsened -> it's all correlated idk
histogram2d(abs.(Jpgrad_at_rs_n2_nrs),deltaprimes_n2_nrs; xlabel="∇Jp (\$\\textrm{A/m}^{2}\$) at rs of 2/1",ylabel = "Δ' of 2/1")
cor(abs.(Jpgrad_at_rs_n2_nrs),deltaprimes_n2_nrs)
corspearman(abs.(Jpgrad_at_rs_n2_nrs),deltaprimes_n2_nrs)

#N3 Jt bout same
histogram2d(abs.(Jtgrad_at_rs_n3_nrs_n3),deltaprimes_n3_nrs_n3; xlabel="∇Jt (\$\\textrm{A/m}^{2}\$) at rs of 2/1",ylabel = "Δ' of 2/1")
cor(abs.(Jtgrad_at_rs_n3_nrs_n3),deltaprimes_n3_nrs_n3)
corspearman(abs.(Jtgrad_at_rs_n3_nrs_n3),deltaprimes_n3_nrs_n3)
#N3 Jp bout same
histogram2d(abs.(Jpgrad_at_rs_n3_nrs_n3),deltaprimes_n3_nrs_n3; xlabel="∇Jp (\$\\textrm{A/m}^{2}\$) at rs of 2/1",ylabel = "Δ' of 2/1")
cor(abs.(Jpgrad_at_rs_n3_nrs_n3),deltaprimes_n3_nrs_n3)
corspearman(abs.(Jpgrad_at_rs_n3_nrs_n3),deltaprimes_n3_nrs_n3)


###################################################### ########################### ########################### 
########################### ########################### ########################### ########################### 
########################### ########################### ########################### ########################### 
#GENERAL HILL:
########################### Second big thing: being on hill correlates with stability. For points in the center of the hill, hill steepness strong correlation with stability
min_closeness=0.0
sep_from_edge=0.05
#hill_equils, a_to_hill, hill_width, hill_gradient, Jt_double_derivative_in_hill, a_to_hill_normed, hill_inds = rs_near_local_minmax(non_rs_equil,min_closeness,sep_from_edge;well=false, return_wellhill_info=true)

cd("/Users/sbenjamin/Desktop/PHD/ScrewPinchDataBase")
#@save "hill_data_n1to7_6677.jld2" hill_equils a_to_hill hill_width hill_gradient Jt_double_derivative_in_hill a_to_hill_normed hill_inds
@load "hill_data_n1to7_6677.jld2" hill_equils a_to_hill hill_width hill_gradient Jt_double_derivative_in_hill a_to_hill_normed hill_inds



#@save "hill_data_n1to4_667.jld2" hill_equils a_to_hill hill_width hill_gradient Jt_double_derivative_in_hill a_to_hill_normed hill_inds
#@load "hill_data_n1to4_667.jld2" hill_equils a_to_hill hill_width hill_gradient Jt_double_derivative_in_hill a_to_hill_normed hill_inds

deltaprime_n2_hill= [i.Δprimes[1].Δprime for i in hill_equils]

hill_gradient[1:10]
plot_current_profiles(hill_equils[1:10];m_ind=1)

min_normed_closeness = 0.1 #Norm to 1
right_in_hill_inds = filter(i -> a_to_hill_normed[i]<min_normed_closeness, 1:length(a_to_hill_normed))
scatter(hill_gradient[right_in_hill_inds],[i.Δprimes[1].Δprime for i in hill_equils[right_in_hill_inds]])
cor(hill_gradient[right_in_hill_inds],[i.Δprimes[1].Δprime for i in hill_equils[right_in_hill_inds]])
corspearman(float.(hill_gradient[right_in_hill_inds]),[i.Δprimes[1].Δprime for i in hill_equils[right_in_hill_inds]])

min_closeness = 0.008 #Stronger correlations, use this
right_in_hill_inds = filter(i -> a_to_hill[i]<min_closeness, 1:length(a_to_hill))
scatter(hill_gradient[right_in_hill_inds],[i.Δprimes[1].Δprime for i in hill_equils[right_in_hill_inds]]; xlabel="∇J/\$\\textrm{J}_{hill}\$", ylabel="Δ' of 2/1", color=:black, label=false)
cor(hill_gradient[right_in_hill_inds],[i.Δprimes[1].Δprime for i in hill_equils[right_in_hill_inds]])
corspearman(float.(hill_gradient[right_in_hill_inds]),[i.Δprimes[1].Δprime for i in hill_equils[right_in_hill_inds]])

right_in_hill_inds_equils=collect(hill_equils[right_in_hill_inds])
right_in_hill_inds_grads=collect(hill_gradient[right_in_hill_inds])
right_in_hill_widths=collect(hill_width[right_in_hill_inds])
Jt_double_derivatives=collect(Jt_double_derivative_in_hill[right_in_hill_inds])
a_to_hill_=collect(a_to_hill[right_in_hill_inds])
a_to_hill_normed_=collect(a_to_hill_normed[right_in_hill_inds])

sorted_hill_equils = collect(right_in_hill_inds_equils[sortperm(right_in_hill_inds_grads)])
sorted_hill_grads = collect(right_in_hill_inds_grads[sortperm(right_in_hill_inds_grads)])
sorted_widths = collect(right_in_hill_widths[sortperm(right_in_hill_inds_grads)])
sorted_Jtdd= collect(Jt_double_derivatives[sortperm(right_in_hill_inds_grads)])
sorted_a_to_hills= collect(a_to_hill_[sortperm(right_in_hill_inds_grads)])
sorted_a_to_hill_normed= collect(a_to_hill_normed_[sortperm(right_in_hill_inds_grads)])

i = (length(sorted_hill_equils)-5):length(sorted_hill_equils)


plot_current_profiles(sorted_hill_equils[end-2];m_ind=nothing, thickness_scaling=1.8,grid=false)#, grid=false,ylims=(2.0,3.0),xlims=(0.4,0.7))

plot_current_profiles(sorted_hill_equils[end-2];m_ind=1, thickness_scaling=1.4, grid=false,ylims=(1.8,3.4),xlims=(0.3,0.8))
plot_current_profiles(sorted_hill_equils[end-2];m_ind=1, thickness_scaling=1.4, grid=false,ylims=(1.8,3.4),xlims=(0.3,0.8))

print("\n grad= $(sorted_hill_grads[i])\n width= $(sorted_widths[i])\n Jt double deriv= $(sorted_Jtdd[i])\n a to min/max= $(sorted_a_to_hills[i])\n a to min/max normed= $(sorted_a_to_hill_normed[i])\n")

    ##########CHOSING M3DC1 runs
        randidxs = rand(1:length(right_in_hill_inds),5)
        scatter(hill_gradient[right_in_hill_inds],[i.Δprimes[1].Δprime for i in hill_equils[right_in_hill_inds]])
        #scatter!(hill_gradient[right_in_hill_inds[randidxs]],[i.Δprimes[1].Δprime for i in hill_equils[right_in_hill_inds[randidxs]]])
        is=[30,49,51,46,79]
        scatter!(hill_gradient[right_in_hill_inds[is]],[i.Δprimes[1].Δprime for i in hill_equils[right_in_hill_inds[is]]])

        #49##
        #51##
        #46##
        #79##
        #30##
        hill_m3dc1_equils = hill_equils[right_in_hill_inds[is]]
        plot_current_profiles(hill_m3dc1_equils[4:5];m_ind=1)
        deltaprime_n2_hill= [i.Δprimes[1].Δprime for i in hill_m3dc1_equils]

########################### What about being in a well? The strongest correlation... why? because well = bad but low current gradient = good, so goes from robustly stable to very unstable as well depth increases from 0. 
#well_equils, a_to_well, well_width, well_gradient, Jt_double_derivative_in_well, a_to_well_normed, well_inds = rs_near_local_minmax(non_rs_equil,min_closeness,sep_from_edge;well=true,return_wellhill_info=true)

###################################################### ########################### ########################### 
########################### ########################### ########################### ########################### 
########################### ########################### ########################### ########################### 
#PLOTTING M3DC1 Cases########################### 
cd("/Users/sbenjamin/Desktop/PHD/ScrewPinchDataBase")
#@save "well_data_n1to7_6677.jld2" well_equils a_to_well well_width well_gradient Jt_double_derivative_in_well a_to_well_normed well_inds
#@load "well_data_n1to7_6677.jld2" well_equils a_to_well well_width well_gradient Jt_double_derivative_in_well a_to_well_normed well_inds

#@save "well_data_n1to4_667.jld2" well_equils a_to_well well_width well_gradient Jt_double_derivative_in_well a_to_well_normed well_inds
@load "well_data_n1to4_667.jld2" well_equils a_to_well well_width well_gradient Jt_double_derivative_in_well a_to_well_normed well_inds

plot_current_profiles(well_equils[1:100];m_ind=1)

deltaprime_n2_well= [i.Δprimes[1].Δprime for i in well_equils]

well_gradient[1:10]
plot_current_profiles(well_equils[1:10];m_ind=1)

min_normed_closeness = 0.1 #Norm to 1
right_in_well_inds = filter(i -> a_to_well_normed[i]<min_normed_closeness, 1:length(a_to_well_normed))
scatter(well_gradient[right_in_well_inds],[i.Δprimes[1].Δprime for i in well_equils[right_in_well_inds]])
cor(well_gradient[right_in_well_inds],[i.Δprimes[1].Δprime for i in well_equils[right_in_well_inds]])
corspearman(float.(well_gradient[right_in_well_inds]),[i.Δprimes[1].Δprime for i in well_equils[right_in_well_inds]])

min_closeness = 0.008 #Stronger correlations, use this
right_in_well_inds = filter(i -> a_to_well[i]<min_closeness, 1:length(a_to_well))
scatter(well_gradient[right_in_well_inds],[i.Δprimes[1].Δprime for i in well_equils[right_in_well_inds]]; xlabel="∇J/\$\\textrm{J}_{well}\$", ylabel="Δ' of 2/1", color=:black, label=false)
cor(well_gradient[right_in_well_inds],[i.Δprimes[1].Δprime for i in well_equils[right_in_well_inds]])
corspearman(float.(well_gradient[right_in_well_inds]),[i.Δprimes[1].Δprime for i in well_equils[right_in_well_inds]])


right_in_well_equils = collect(well_equils[right_in_well_inds])
right_in_well_gradient = collect(well_gradient[right_in_well_inds])

sorted_well_equils667 = collect(right_in_well_equils[sortperm([i.Δprimes[1].Δprime for i in right_in_well_equils])])
sorted_well_grad667 = collect(right_in_well_gradient[sortperm([i.Δprimes[1].Δprime for i in right_in_well_equils])])

i=4
scatter(sorted_well_grad667,[i.Δprimes[1].Δprime for i in sorted_well_equils667])
scatter!(sorted_well_grad667[i:i],[j.Δprimes[1].Δprime for j in sorted_well_equils667[i:i]]; label=false)

iswell=[4,41,60,80,100,150,230,284]
well_m3dc1_equils=collect(sorted_well_equils667[iswell])
well_m3dc1_grads=collect(sorted_well_grad667[iswell])
plot_current_profiles(well_m3dc1_equils;m_ind=1)
deltaprime_n2_well= [i.Δprimes[1].Δprime for i in well_m3dc1_equils]
well_m3dc1_equils_extra=sorted_well_equils667[1]
iswellextrs=[1,4,41,60,80,100,150,230,284]
Jtconfirms = [j.equilibrium.Jt(0.01) for j in sorted_well_equils667[iswellextrs]]

scatter(sorted_well_grad667,[i.Δprimes[1].Δprime for i in sorted_well_equils667];label=false)
scatter!(sorted_well_grad667[iswellextrs],[j.Δprimes[1].Δprime for j in sorted_well_equils667[iswellextrs]]; label=false)

plot_current_profiles(sorted_well_equils667[iswellextrs[1:3]];m_ind=1)


sorted_well_dpthequils = collect(right_in_well_equils[sortperm(right_in_well_gradient)])
sorted_well_dpthgrad = collect(sort(right_in_well_gradient))

plot_current_profiles(sorted_well_dpthequils[end-4];m_ind=nothing, thickness_scaling=1.8,grid=false)#, grid=false,ylims=(2.0,3.0),xlims=(0.4,0.7))

plot_current_profiles(sorted_well_dpthequils[end-4];m_ind=1, thickness_scaling=1.4, grid=false,ylims=(0,3),xlims=(0.3,0.95))


###################################################### ########################### ########################### 
########################### ########################### ########################### ########################### 
########################### ########################### ########################### ########################### 
#GENERAL WELL
cd("/Users/sbenjamin/Desktop/PHD/ScrewPinchDataBase")
#@save "well_data_n1to7_6677.jld2" well_equils a_to_well well_width well_gradient Jt_double_derivative_in_well a_to_well_normed well_inds
@load "well_data_n1to7_6677.jld2" well_equils a_to_well well_width well_gradient Jt_double_derivative_in_well a_to_well_normed well_inds

plot_current_profiles(well_equils[1:100];m_ind=1)

deltaprime_n2_well= [i.Δprimes[1].Δprime for i in well_equils]

well_gradient[1:10]
plot_current_profiles(well_equils[1:10];m_ind=1)

min_normed_closeness = 0.1 #Norm to 1
right_in_well_inds = filter(i -> a_to_well_normed[i]<min_normed_closeness, 1:length(a_to_well_normed))
scatter(well_gradient[right_in_well_inds],[i.Δprimes[1].Δprime for i in well_equils[right_in_well_inds]])
cor(well_gradient[right_in_well_inds],[i.Δprimes[1].Δprime for i in well_equils[right_in_well_inds]])
corspearman(float.(well_gradient[right_in_well_inds]),[i.Δprimes[1].Δprime for i in well_equils[right_in_well_inds]])

min_closeness = 0.008 #Stronger correlations, use this
right_in_well_inds = filter(i -> a_to_well[i]<min_closeness, 1:length(a_to_well))
scatter(well_gradient[right_in_well_inds],[i.Δprimes[1].Δprime for i in well_equils[right_in_well_inds]]; xlabel="∇J/\$\\textrm{J}_{well}\$", ylabel="Δ' of 2/1", color=:black, label=false)
cor(well_gradient[right_in_well_inds],[i.Δprimes[1].Δprime for i in well_equils[right_in_well_inds]])
corspearman(float.(well_gradient[right_in_well_inds]),[i.Δprimes[1].Δprime for i in well_equils[right_in_well_inds]])
scatter!(well_m3dc1_grads[2:end],[j.Δprimes[1].Δprime for j in well_m3dc1_equils[2:end]]; xlabel="∇J/\$\\textrm{J}_{well}\$", ylabel="Δ' of 2/1", color=:red, label=false)

"""ååå
    right_in_well_equils = collect(well_equils[right_in_well_inds])
    right_in_well_gradient = collect(well_gradient[right_in_well_inds])

    sorted_well_equils = collect(right_in_well_equils[sortperm([i.Δprimes[1].Δprime for i in right_in_well_equils])])
    sorted_well_grad = collect(right_in_well_gradient[sortperm([i.Δprimes[1].Δprime for i in right_in_well_equils])])

    i=4
    scatter(sorted_well_grad,[i.Δprimes[1].Δprime for i in sorted_well_equils])
    scatter!(sorted_well_grad[i:i],[j.Δprimes[1].Δprime for j in sorted_well_equils[i:i]]; label=false)

    iswell=[4,41,60,80,100,150,230,284]
    well_m3dc1_equils=collect(sorted_well_equils[iswell])
    plot_current_profiles(well_m3dc1_equils;m_ind=1)
    deltaprime_n2_well= [i.Δprimes[1].Δprime for i in well_m3dc1_equils]
    well_m3dc1_equils_extra=sorted_well_equils[1]
    iswellextrs=[1,4,41,60,80,100,150,230,284]
"""

scatter(sorted_well_grad,[i.Δprimes[1].Δprime for i in sorted_well_equils];label=false)
scatter!(sorted_well_grad667[iswellextrs],[j.Δprimes[1].Δprime for j in sorted_well_equils667[iswellextrs]]; label=false)

Jtconfirms = [j.equilibrium.Jt(0.0) for j in sorted_well_equils667[iswellextrs]]


right_in_well_equils = collect(well_equils[right_in_well_inds])
right_in_well_gradient = collect(well_gradient[right_in_well_inds])

###################################################### ########################### ########################### 
########################### ########################### ########################### ########################### 
########################### ########################### ########################### ########################### 
#OUT WELL:
welll_equils, a_to_welll, welll_width, welll_gradient, Jt_double_derivative_in_welll, a_to_welll_normed, welll_inds, welll_edge_vals = rs_near_local_minmax(right_in_well_equils,min_closeness,sep_from_edge;well=true,return_wellhill_info=true, extra_outputs=true)

∇J_out_wrtrs = []
∇J_out_wrtwell = []
for i in welll_edge_vals
    push!(∇J_out_wrtrs,((1/(100*100))*(i[1,2]-i[3,2])/(i[1,1]-i[3,1]))/((1/(100*100))*i[3,2]))
    push!(∇J_out_wrtwell,((1/(100*100))*(i[1,2]-i[2,2])/(i[1,1]-i[2,1]))/((1/(100*100))*i[2,2]))
    if i[1,2]-i[3,2]<0 || i[1,1]-i[3,1]<0 || i[1,2]-i[2,2]<0 || i[1,1]-i[2,1]<0
        print("wth")
    end
end

scatter(∇J_out_wrtrs,[i.Δprimes[1].Δprime for i in welll_equils])
scatter(∇J_out_wrtwell,[i.Δprimes[1].Δprime for i in welll_equils]; xlabel="∇J/\$\\textrm{J}_{out}\$", ylabel="Δ' of 2/1", color=:black, label=false, thickness_scaling=1.5)
corspearman(float.(∇J_out_wrtrs),[i.Δprimes[1].Δprime for i in welll_equils])
corspearman(float.(∇J_out_wrtwell),[i.Δprimes[1].Δprime for i in welll_equils])





########################### Some general correlations:
Jts_nrs = [i.equilibrium.Jt for i in non_rs_equil]                                          #
Jps_nrs = [i.equilibrium.Jp for i in non_rs_equil]                                          #
rss_n2_nrs = [i.Δprimes[1].rs for i in non_rs_equil]                                        #
rss_n3_nrs = [i.Δprimes[2].rs for i in non_rs_equil]                                        #
Jtots_nrs = [total_plasma_current(i.equilibrium.Jt,i.equilibrium.rb) for i in non_rs_equil] #
Jt_at_rs_n2_nrs = [Jts_nrs[i](rss_n2_nrs[i]) for i in 1:length(non_rs_equil)]               
Jt_at_rs_n3_nrs = [Jts_nrs[i](rss_n3_nrs[i]) for i in 1:length(non_rs_equil)]
Jp_at_rs_n2_nrs = [Jps_nrs[i](rss_n2_nrs[i]) for i in 1:length(non_rs_equil)]
Jp_at_rs_n3_nrs = [Jps_nrs[i](rss_n3_nrs[i]) for i in 1:length(non_rs_equil)]
shear_at_rs_n2_nrs = [ForwardDiff.derivative(non_rs_equil[i].equilibrium.q,rss_n2_nrs[i]) for i in 1:length(non_rs_equil)]
shear_at_rs_n3_nrs = [ForwardDiff.derivative(non_rs_equil[i].equilibrium.q,rss_n3_nrs[i]) for i in 1:length(non_rs_equil)]
qdd_at_rs_n2_nrs = [ForwardDiff.derivative(x0->ForwardDiff.derivative(non_rs_equil[i].equilibrium.q,x0),rss_n2_nrs[i]) for i in 1:length(non_rs_equil)]
q_at_rs_n2_nrs = [non_rs_equil[i].equilibrium.q(rss_n2_nrs[i]) for i in 1:length(non_rs_equil)]
q_at_rs_n3_nrs = [non_rs_equil[i].equilibrium.q(rss_n3_nrs[i]) for i in 1:length(non_rs_equil)]
p_at_rs_n2_nrs= [mu0*non_rs_equil[i].equilibrium.p(rss_n2_nrs[i]) for i in 1:length(non_rs_equil)]
p_at_rs_n3_nrs = [mu0*non_rs_equil[i].equilibrium.p(rss_n3_nrs[i]) for i in 1:length(non_rs_equil)]
dpdr_at_rs_n2_nrs = [non_rs_equil[i].equilibrium.dpdr(rss_n2_nrs[i]) for i in 1:length(non_rs_equil)]
dpdr_at_rs_n3_nrs = [non_rs_equil[i].equilibrium.dpdr(rss_n3_nrs[i]) for i in 1:length(non_rs_equil)]



histogram2d(Jtots_nrs./1e6,deltaprimes_n2_nrs; thickness_scaling=1.4, xlabel="Total current (MA)",ylabel = "Δ' 2/1") #NO CORRELATION
cor(Jtots_nrs,deltaprimes_n2_nrs)
corspearman(Jtots_nrs,deltaprimes_n2_nrs)
histogram2d(Jt_at_rs_n2_nrs./1e6,deltaprimes_n2_nrs;thickness_scaling=1.4, xlabel="Jt at rs 2/1 (\$\\textrm{MA/m}^{2}\$)",ylabel = "Δ' 2/1") #SMALL CORRELATION
cor(Jt_at_rs_n2_nrs,deltaprimes_n2_nrs)
corspearman(Jt_at_rs_n2_nrs,deltaprimes_n2_nrs)
histogram2d(Jp_at_rs_n2_nrs./1e6,deltaprimes_n2_nrs;thickness_scaling=1.4, xlabel="Jp at rs 2/1 (\$\\textrm{MA/m}^{2}\$)",ylabel = "Δ' 2/1") #ALMOST NO CORRELATION
cor(Jp_at_rs_n2_nrs,deltaprimes_n2_nrs)
corspearman(Jp_at_rs_n2_nrs,deltaprimes_n2_nrs) 
histogram2d(p_at_rs_n2_nrs,deltaprimes_n2_nrs;thickness_scaling=1.3, xlabel="Pressure at rs 2/1 (\$\\sqrt{T}\$)",ylabel = "Δ' 2/1")  #SMALL CORRELATION (parabolic pressure profile its about location of rs)
cor(p_at_rs_n2_nrs,deltaprimes_n2_nrs) 
corspearman(p_at_rs_n2_nrs,deltaprimes_n2_nrs) 
histogram2d(dpdr_at_rs_n2_nrs,deltaprimes_n2_nrs; thickness_scaling=1.3,xlabel="Pressure gradient at rs 2/1 (\$\\sqrt{T}\$/m)",ylabel = "Δ' 2/1")   #SMALL CORRELATION LIKE p (parabolic pressure profile its about location of rs)
cor(dpdr_at_rs_n2_nrs,deltaprimes_n2_nrs)
corspearman(dpdr_at_rs_n2_nrs,deltaprimes_n2_nrs)
histogram2d(rss_n2_nrs,deltaprimes_n2_nrs; thickness_scaling=1.3,xlabel="rs 2/1 (sqrt(T)/m)",ylabel = "Δ' 2/1")   #SMALL CORRELATION LIKE p (parabolic pressure profile its about location of rs)
cor(rss_n2_nrs,deltaprimes_n2_nrs)
corspearman(rss_n2_nrs,deltaprimes_n2_nrs)

histogram2d(Jpgrad_at_rs_n2_nrs./1e6,deltaprimes_n2_nrs;  thickness_scaling=1.4, xlabel="Jp' (\$\\textrm{MA/m}^{3}\$) at rs of 2/1",ylabel = "Δ' of 2/1")
cor(Jpgrad_at_rs_n2_nrs,deltaprimes_n2_nrs)
corspearman(Jpgrad_at_rs_n2_nrs,deltaprimes_n2_nrs)


histogram2d(q_at_rs_n2_nrs,deltaprimes_n2_nrs; xlabel="q at rs 2/1",ylabel = "Δ' 2/1") #SMALL CORRELATION
cor(q_at_rs_n2_nrs,deltaprimes_n2_nrs)
corspearman(q_at_rs_n2_nrs,deltaprimes_n2_nrs)

histogram2d(shear_at_rs_n2_nrs,deltaprimes_n2_nrs; thickness_scaling=1.4, xlabel="q' at rs 2/1",ylabel = "Δ' 2/1") #SMALL CORRELATION
cor(shear_at_rs_n2_nrs,deltaprimes_n2_nrs)
corspearman(shear_at_rs_n2_nrs,deltaprimes_n2_nrs)

histogram2d(qdd_at_rs_n2_nrs,deltaprimes_n2_nrs; thickness_scaling=1.4, xlabel="q'' at rs 2/1",ylabel = "Δ' 2/1") #SMALL CORRELATION
cor(qdd_at_rs_n2_nrs,deltaprimes_n2_nrs)
corspearman(qdd_at_rs_n2_nrs,deltaprimes_n2_nrs)

histogram2d(Jtgrad_at_rs_n2_nrs./1e6,deltaprimes_n2_nrs; thickness_scaling=1.25, xlabel="Jt' (\$\\textrm{MA/m}^{3}\$) at rs of 2/1",ylabel = "Δ' of 2/1")
cor(Jtgrad_at_rs_n2_nrs,deltaprimes_n2_nrs)
corspearman(Jtgrad_at_rs_n2_nrs,deltaprimes_n2_nrs)

histogram2d(Jt_at_rs_n2_nrs./1e6,deltaprimes_n2_nrs;thickness_scaling=1.4, xlabel="Jt at rs 2/1 (\$\\textrm{MA/m}^{2}\$)",ylabel = "Δ' 2/1") #SMALL CORRELATION
cor(Jt_at_rs_n2_nrs,deltaprimes_n2_nrs)
corspearman(Jt_at_rs_n2_nrs,deltaprimes_n2_nrs)


############################## REMOVING Wells and Hills:
########################### Some general correlations:

non_rs_equil_inds_no_wells = Int[]
non_rs_equil_inds_no_hills = Int[]
non_rs_equil_inds_no_well_hills = Int[]

for i in 1:length(non_rs_equil)
    if !(i in well_inds)
        push!(non_rs_equil_inds_no_wells,i)
    end
    if !(i in hill_inds)
        push!(non_rs_equil_inds_no_hills,i)
    end
    if !(i in hill_inds) && !(i in well_inds)
        push!(non_rs_equil_inds_no_well_hills,i)
    end
end

non_rs_equil_nowells = collect(non_rs_equil[non_rs_equil_inds_no_wells])
non_rs_equil_nowellhills = collect(non_rs_equil[non_rs_equil_inds_no_well_hills])

Jts_nrs_nowells = [i.equilibrium.Jt for i in non_rs_equil_nowells] 
Jts_nrs_nowellhills = [i.equilibrium.Jt for i in non_rs_equil_nowellhills] 

rss_n2_nrs_nowells = [i.Δprimes[1].rs for i in non_rs_equil_nowells]   
rss_n2_nrs_nowellhills = [i.Δprimes[1].rs for i in non_rs_equil_nowellhills]   

Jtgrad_at_rs_n2_nrs_nowells = [ScrewPinchAsymptoticMatching.gradient(Jts_nrs_nowells[i],rss_n2_nrs_nowells[i]) for i in 1:length(non_rs_equil_nowells)]
Jtgrad_at_rs_n2_nrs_nowellhills = [ScrewPinchAsymptoticMatching.gradient(Jts_nrs_nowellhills[i],rss_n2_nrs_nowellhills[i]) for i in 1:length(non_rs_equil_nowellhills)]
Jt_at_rs_n2_nrs_nowellhills = [Jts_nrs_nowellhills[i](rss_n2_nrs_nowellhills[i]) for i in 1:length(non_rs_equil_nowellhills)]

deltaprimes_n2_nrs_nowells= [i.Δprimes[1].Δprime for i in non_rs_equil_nowells]
deltaprimes_n2_nrs_nowellhills= [i.Δprimes[1].Δprime for i in non_rs_equil_nowellhills]

sorted_inds_nowell = sortperm(deltaprimes_n2_nrs_nowells)
sorted_inds_nowellhill = sortperm(deltaprimes_n2_nrs_nowellhills)

deltaprimes_n2_nrs_nowells[sorted_inds_nowell[1:10]]
deltaprimes_n2_nrs_nowellhills[sorted_inds_nowellhill[1:10]]

deltaprimes_n2_nrs_nowells[sorted_inds_nowell[end-10:end]]
deltaprimes_n2_nrs_nowellhills[sorted_inds_nowellhill[end-10:end]]

#orig
    #Jt n2 big correlation
    histogram2d(Jtgrad_at_rs_n2_nrs./1e6,deltaprimes_n2_nrs; xlabel="\$J_t^'(r_s) (\\text{MA/m}^{3}) \\textrm{at } \\textrm{r_s}\$ of 2/1",ylabel = "Δ' of 2/1")
    cor(Jtgrad_at_rs_n2_nrs,deltaprimes_n2_nrs)
    corspearman(Jtgrad_at_rs_n2_nrs,deltaprimes_n2_nrs)

#nowells
    #Jt n2 big correlation
    histogram2d(Jtgrad_at_rs_n2_nrs_nowells./1e6,deltaprimes_n2_nrs_nowells; xlabel="\$J_t^'(r_s) (\\text{MA/m}^{3}) \\textrm{at } \\textrm{r_s}\$ of 2/1",ylabel = "Δ' of 2/1")
    cor(Jtgrad_at_rs_n2_nrs_nowells,deltaprimes_n2_nrs_nowells)
    corspearman(Jtgrad_at_rs_n2_nrs_nowells,deltaprimes_n2_nrs_nowells)

#nowellhills
    #Jt n2 big correlation
    histogram2d(Jtgrad_at_rs_n2_nrs_nowellhills./1e6,deltaprimes_n2_nrs_nowellhills; xlabel="\$J_t^'(r_s) (\\text{MA/m}^{3}) \\textrm{at } \\textrm{r_s}\$ of 2/1",ylabel = "Δ' of 2/1")
    cor(Jtgrad_at_rs_n2_nrs_nowellhills,deltaprimes_n2_nrs_nowellhills)
    corspearman(Jtgrad_at_rs_n2_nrs_nowellhills,deltaprimes_n2_nrs_nowellhills)

############################## Final figs
########################### 

cdog=:tofino
clin=cgrad(cdog, scale = :linear)
clog=cgrad(cdog, scale = :exp)


histogram2d(rss_n2,deltaprimes_n2; c=clog,xlabel="2/1 rational surface location (m)",ylabel = "Δ' 2/1")
cor(rss_n2,deltaprimes_n2)
corspearman(rss_n2,deltaprimes_n2)
savefig()
histogram2d(rss_n3[n_3_inds],deltaprimes_n3[n_3_inds]; c=clog,xlabel="3/1 rational surface location (m)",ylabel = "Δ' 3/1") 
cor(rss_n3[n_3_inds],deltaprimes_n3[n_3_inds])
corspearman(rss_n3[n_3_inds],deltaprimes_n3[n_3_inds])
savefig[]


histogram2d(Jtgrad_at_rs_n2_nrs./1e6,deltaprimes_n2_nrs;  c=clog, xlabel="\$J_t^'(r_s) (\\text{MA/m}^{3}) \\textrm{at } \\textrm{r_s}\$ of 2/1",ylabel = "Δ' of 2/1")
cor(Jtgrad_at_rs_n2_nrs,deltaprimes_n2_nrs)
corspearman(Jtgrad_at_rs_n2_nrs,deltaprimes_n2_nrs)
savefig[]

histogram2d(Jtgrad_at_rs_n2_nrs./1e6,deltaprimes_n2_nrs; c=clog, xrange=[-Inf,0.3+maximum(Jtgrad_at_rs_n2_nrs)/1e6],xlabel="\$J_t^'(r_s) (\\text{MA/m}^{3}) \\textrm{at } \\textrm{r_s}\$ of 2/1",ylabel = "Δ' of 2/1")#, xlabel="\$J_t'(r_s)\$ (MA/m\$^3\$) at r\$_s\$ of \$2/1\$",ylabel = "Δ\$'\$ of 2/1",fontfamily="serif-roman")
cor(Jtgrad_at_rs_n2_nrs,deltaprimes_n2_nrs)
corspearman(Jtgrad_at_rs_n2_nrs,deltaprimes_n2_nrs)
savefig[]
histogram2d(Jtgrad_at_rs_n2_nrs_nowells./1e6,deltaprimes_n2_nrs_nowells; c=clog, xrange=[-Inf,0.3+maximum(Jtgrad_at_rs_n2_nrs)/1e6], xlabel="\$J_t^'(r_s) (\\text{MA/m}^{3}) \\textrm{at } \\textrm{r_s}\$ of 2/1",ylabel = "Δ' of 2/1")
cor(Jtgrad_at_rs_n2_nrs_nowells,deltaprimes_n2_nrs_nowells)
corspearman(Jtgrad_at_rs_n2_nrs_nowells,deltaprimes_n2_nrs_nowells)
savefig[]
histogram2d(Jtgrad_at_rs_n2_nrs_nowellhills./1e6,deltaprimes_n2_nrs_nowellhills; c=clog, yrange=[-27.5,Inf], xrange=[-Inf,0.3+maximum(Jtgrad_at_rs_n2_nrs)/1e6], xlabel="\$J_t^'(r_s) (\\text{MA/m}^{3}) \\textrm{at } \\textrm{r_s}\$ of 2/1",ylabel = "Δ' of 2/1")
cor(Jtgrad_at_rs_n2_nrs_nowellhills,deltaprimes_n2_nrs_nowellhills)
corspearman(Jtgrad_at_rs_n2_nrs_nowellhills,deltaprimes_n2_nrs_nowellhills)
savefig
histogram2d(filter(x->(x>-8)&&(x<0.2),collect(Jtgrad_at_rs_n2_nrs_nowells./1e6)),deltaprimes_n2_nrs_nowells[filter(x->Jtgrad_at_rs_n2_nrs_nowells[x]/1e6>-8&&Jtgrad_at_rs_n2_nrs_nowells[x]/1e6<0.2,1:length(Jtgrad_at_rs_n2_nrs_nowells))]; bins=(50,500), c=clog, xrange=[-8,0.2],yrange=[-9,19],label=nothing,axes=false,thickness_scaling=1.9,framestyle = :box)
cor(Jtgrad_at_rs_n2_nrs_nowells,deltaprimes_n2_nrs_nowells)
corspearman(Jtgrad_at_rs_n2_nrs_nowells,deltaprimes_n2_nrs_nowells)
xdots=StepRangeLen(-7.5,0.1,60)
plot!(xdots,xdots.*(-2.15).-10,linewidth=8,linecolor=:red,line=:dot,label=false)

histogram2d(filter(x->(x>-8)&&(x<0.2),collect(Jtgrad_at_rs_n2_nrs_nowells./1e6)),deltaprimes_n2_nrs_nowells[filter(x->Jtgrad_at_rs_n2_nrs_nowells[x]/1e6>-8&&Jtgrad_at_rs_n2_nrs_nowells[x]/1e6<0.2,1:length(Jtgrad_at_rs_n2_nrs_nowells))]; bins=(50,500), c=clog, xrange=[-8,0.2],yrange=[-9,19],label=nothing,axes=false,thickness_scaling=1.9,ticks=false,framestyle = :box)
cor(Jtgrad_at_rs_n2_nrs_nowells,deltaprimes_n2_nrs_nowells)
corspearman(Jtgrad_at_rs_n2_nrs_nowells,deltaprimes_n2_nrs_nowells)
#xdots=StepRangeLen(-8,0.1,100)
#plot!(xdots,xdots.*0.0,linewidth=4,linecolor=:black,line=:line,label=false)
xdots=StepRangeLen(-7.5,0.1,60)
plot!(xdots,xdots.*(-2.15).-10,linewidth=8, xrange=[-8,0.2],linecolor=:red,line=:dot,label=false)

histogram2d(filter(x->(x>-8)&&(x<0.2),collect(Jtgrad_at_rs_n2_nrs_nowellhills./1e6)),deltaprimes_n2_nrs_nowellhills[filter(x->Jtgrad_at_rs_n2_nrs_nowellhills[x]/1e6>-8&&Jtgrad_at_rs_n2_nrs_nowellhills[x]/1e6<0.2,1:length(Jtgrad_at_rs_n2_nrs_nowellhills))]; bins=(50,500), c=clog, xrange=[-8,0.2],yrange=[-9,19],label=nothing,axes=false,thickness_scaling=1.9,framestyle = :box)
cor(Jtgrad_at_rs_n2_nrs_nowellhills,deltaprimes_n2_nrs_nowellhills)
corspearman(Jtgrad_at_rs_n2_nrs_nowellhills,deltaprimes_n2_nrs_nowellhills)
#xdots=StepRangeLen(-8,0.1,100)
#plot!(xdots,xdots.*0.0,linewidth=4,linecolor=:black,line=:line,label=false)
xdots=StepRangeLen(-7.5,0.1,60)
plot!(xdots,xdots.*(-2.15).-10,linewidth=8, xrange=[-8,0.2],linecolor=:red,line=:dot,label=false)
xdots2=StepRangeLen(-4.5,0.1,43)
plot!(xdots2,xdots2.*(-7.15).-7,linewidth=8, xrange=[-8,0.2],linecolor=:red,line=:dot,label=false)

histogram2d(filter(x->x<0,collect(Jtgrad_at_rs_n2_nrs_nowellhills./1e6)),deltaprimes_n2_nrs_nowellhills[filter(x->Jtgrad_at_rs_n2_nrs_nowellhills[x]/1e6<0,1:length(Jtgrad_at_rs_n2_nrs_nowellhills))])
cor(filter(x->x<0,collect(Jtgrad_at_rs_n2_nrs_nowellhills./1e6)),deltaprimes_n2_nrs_nowellhills[filter(x->Jtgrad_at_rs_n2_nrs_nowellhills[x]/1e6<0,1:length(Jtgrad_at_rs_n2_nrs_nowellhills))])
corspearman(filter(x->x<0,collect(Jtgrad_at_rs_n2_nrs_nowellhills./1e6)),deltaprimes_n2_nrs_nowellhills[filter(x->Jtgrad_at_rs_n2_nrs_nowellhills[x]/1e6<0,1:length(Jtgrad_at_rs_n2_nrs_nowellhills))])

savefig[]


histogram2d(shear_at_rs_n2_nrs,deltaprimes_n2_nrs; c=clog,thickness_scaling=1.4, xlabel="q' at rs 2/1",ylabel = "Δ' 2/1",size = (1.5*500,600)) #SMALL CORRELATION
cor(shear_at_rs_n2_nrs,deltaprimes_n2_nrs)
corspearman(shear_at_rs_n2_nrs,deltaprimes_n2_nrs)
savefig[]
histogram2d(qdd_at_rs_n2_nrs,deltaprimes_n2_nrs; c=clog,thickness_scaling=1.4, xlabel="q'' at rs 2/1",ylabel = "Δ' 2/1",size = (1.5*500,600))#, right_margin =7Plots.px) #SMALL CORRELATION
cor(qdd_at_rs_n2_nrs,deltaprimes_n2_nrs)
corspearman(qdd_at_rs_n2_nrs,deltaprimes_n2_nrs)
savefig[]
histogram2d(Jtgrad_at_rs_n2_nrs./1e6,deltaprimes_n2_nrs; c=clog,thickness_scaling=1.25, xlabel="Jt' (\$\\textrm{MA/m}^{3}\$) at rs of 2/1",ylabel = "Δ' of 2/1",size = (1.5*500,600))
cor(Jtgrad_at_rs_n2_nrs,deltaprimes_n2_nrs)
corspearman(Jtgrad_at_rs_n2_nrs,deltaprimes_n2_nrs)
savefig[]
histogram2d(Jt_at_rs_n2_nrs./1e6,deltaprimes_n2_nrs;c=clog,thickness_scaling=1.4, xlabel="Jt at rs 2/1 (\$\\textrm{MA/m}^{2}\$)",ylabel = "Δ' 2/1",size = (1.5*500,600)) #SMALL CORRELATION
cor(Jt_at_rs_n2_nrs,deltaprimes_n2_nrs)
corspearman(Jt_at_rs_n2_nrs,deltaprimes_n2_nrs)
savefig[]


#NORMED SOLNS w wells, hills
histogram2d(Jtgrad_at_rs_n2_nrs./Jt_at_rs_n2_nrs,deltaprimes_n2_nrs; c=clog,thickness_scaling=1.25, xlabel="Jt' (\$\\textrm{MA/m}^{3}\$) at rs of 2/1",ylabel = "Δ' of 2/1",size = (1.5*500,600))
cor(Jtgrad_at_rs_n2_nrs./Jt_at_rs_n2_nrs,deltaprimes_n2_nrs)
corspearman(Jtgrad_at_rs_n2_nrs./Jt_at_rs_n2_nrs,deltaprimes_n2_nrs)
#NORMED SOLNS w no wells, hills, no rs-scaling
histogram2d(Jtgrad_at_rs_n2_nrs_nowellhills./Jt_at_rs_n2_nrs_nowellhills,deltaprimes_n2_nrs_nowellhills; c=clog,thickness_scaling=1.25, xlabel="Jt'/Jt (\$\\textrm{MA/m}^{3}\$) at rs of 2/1",ylabel = "Δ' of 2/1",size = (1.5*500,600))
cor(Jtgrad_at_rs_n2_nrs_nowellhills./Jt_at_rs_n2_nrs_nowellhills,deltaprimes_n2_nrs_nowellhills)
corspearman(Jtgrad_at_rs_n2_nrs_nowellhills./Jt_at_rs_n2_nrs_nowellhills,deltaprimes_n2_nrs_nowellhills)

#NORMED SOLNS w no wells, hills
#rs*Jt'/Jt vs rs*Delta'
histogram2d(rss_n2_nrs_nowellhills.*(Jtgrad_at_rs_n2_nrs_nowellhills./Jt_at_rs_n2_nrs_nowellhills),rss_n2_nrs_nowellhills.*deltaprimes_n2_nrs_nowellhills; c=clog,thickness_scaling=1.25, bins=(200,200), xlabel="rs(Jt'/Jt) (\$\\textrm{MA/m}^{3}\$) at rs of 2/1",ylabel = "rsΔ' of 2/1",size = (1.5*500,600))
cor(rss_n2_nrs_nowellhills.*(Jtgrad_at_rs_n2_nrs_nowellhills./Jt_at_rs_n2_nrs_nowellhills),rss_n2_nrs_nowellhills.*deltaprimes_n2_nrs_nowellhills)
corspearman(rss_n2_nrs_nowellhills.*(Jtgrad_at_rs_n2_nrs_nowellhills./Jt_at_rs_n2_nrs_nowellhills),rss_n2_nrs_nowellhills.*deltaprimes_n2_nrs_nowellhills)



goofynormedJtgrads = collect(rss_n2_nrs_nowellhills.*(Jtgrad_at_rs_n2_nrs_nowellhills./Jt_at_rs_n2_nrs_nowellhills))
goofynormeddeltys = collect(rss_n2_nrs_nowellhills.*deltaprimes_n2_nrs_nowellhills)
#normedJtgrads = filter(i->i>-4, goofynormedJtgrads)
normedJtgrads = collect(goofynormedJtgrads[filter(i->goofynormedJtgrads[i]>-3, 1:length(goofynormedJtgrads))])
normeddeltys = collect(goofynormeddeltys[filter(i->goofynormedJtgrads[i]>-3, 1:length(goofynormedJtgrads))])
#normeddeltys = collect(goofynormeddeltys[normeddetlysind])

#NORMED SOLNS w no wells, hills
histogram2d(rss_n2_nrs_nowellhills.*(Jtgrad_at_rs_n2_nrs_nowellhills./Jt_at_rs_n2_nrs_nowellhills),rss_n2_nrs_nowellhills.*deltaprimes_n2_nrs_nowellhills; c=clog,thickness_scaling=1.25, bins=(200,200), xlabel="rs(Jt'/Jt) (\$\\textrm{MA/m}^{3}\$) at rs of 2/1",ylabel = "rsΔ' of 2/1",size = (1.5*500,600))
cor(rss_n2_nrs_nowellhills.*(Jtgrad_at_rs_n2_nrs_nowellhills./Jt_at_rs_n2_nrs_nowellhills),rss_n2_nrs_nowellhills.*deltaprimes_n2_nrs_nowellhills)
corspearman(rss_n2_nrs_nowellhills.*(Jtgrad_at_rs_n2_nrs_nowellhills./Jt_at_rs_n2_nrs_nowellhills),rss_n2_nrs_nowellhills.*deltaprimes_n2_nrs_nowellhills)
#NORMED SOLNS w no wells, hills, DETAIL
histogram2d(normedJtgrads,normeddeltys;range=(-3,6),xlims=(-3,0.1), c=clog,thickness_scaling=1.25, bins=(500,1000), xlabel="rs(Jt'/Jt) (\$\\textrm{MA/m}^{3}\$) at rs of 2/1",ylabel = "rsΔ' of 2/1",size = (1.5*500,600))
cor(normedJtgrads,normeddeltys)
corspearman(normedJtgrads,normeddeltys)


histogram2d(p_at_rs_n2_nrs,deltaprimes_n2_nrs;dpi=1200,c=clog,thickness_scaling=1.8, xlabel="Pressure at rs 2/1 (\$\\sqrt{T}\$)",ylabel = "Δ' 2/1",size = (2*500,2*400))  #SMALL CORRELATION (parabolic pressure profile its about location of rs)
cor(p_at_rs_n2_nrs,deltaprimes_n2_nrs) 
corspearman(p_at_rs_n2_nrs,deltaprimes_n2_nrs) 
savefig("/Users/sbenjamin/Desktop/p.png")
histogram2d(dpdr_at_rs_n2_nrs,deltaprimes_n2_nrs;dpi=1200,c=clog, thickness_scaling=1.8,xlabel="Pressure gradient at rs 2/1 (\$\\sqrt{T}\$/m)",ylabel = "Δ' 2/1",size = (2*500,2*400))    #SMALL CORRELATION LIKE p (parabolic pressure profile its about location of rs)
cor(dpdr_at_rs_n2_nrs,deltaprimes_n2_nrs)
corspearman(dpdr_at_rs_n2_nrs,deltaprimes_n2_nrs)
savefig("/Users/sbenjamin/Desktop/pprime.png")
histogram2d(rss_n2_nrs,deltaprimes_n2_nrs;dpi=1200, c=clog,thickness_scaling=1.8,xlabel="rs 2/1 (sqrt(T)/m)",ylabel = "Δ' 2/1",size = (2*500,2*400))    #SMALL CORRELATION LIKE p (parabolic pressure profile its about location of rs)
cor(rss_n2_nrs,deltaprimes_n2_nrs)
corspearman(rss_n2_nrs,deltaprimes_n2_nrs)
savefig("/Users/sbenjamin/Desktop/rs.png")

histogram2d(Jp_at_rs_n2_nrs./1e6,deltaprimes_n2_nrs;dpi=1200,c=clog,thickness_scaling=1.8, xlabel="Jp at rs 2/1 (\$\\textrm{MA/m}^{2}\$)",ylabel = "Δ' 2/1",size = (2*500,2*400))   #ALMOST NO CORRELATION
cor(Jp_at_rs_n2_nrs,deltaprimes_n2_nrs)
corspearman(Jp_at_rs_n2_nrs,deltaprimes_n2_nrs) 
savefig("/Users/sbenjamin/Desktop/Jp.png")
histogram2d(Jpgrad_at_rs_n2_nrs./1e6,deltaprimes_n2_nrs;dpi=1200,c=clog,  thickness_scaling=1.8, xlabel="Jp' (\$\\textrm{MA/m}^{3}\$) at rs of 2/1",ylabel = "Δ' of 2/1",size = (2*500,2*400)) 
cor(Jpgrad_at_rs_n2_nrs,deltaprimes_n2_nrs)
corspearman(Jpgrad_at_rs_n2_nrs,deltaprimes_n2_nrs)
savefig("/Users/sbenjamin/Desktop/Jpprime.png")
histogram2d(Jtots_nrs./1e6,deltaprimes_n2_nrs;dpi=1200,c=clog, thickness_scaling=1.8, xlabel="Total current (MA)",ylabel = "Δ' 2/1",bins=(200,200),size = (2*500,2*400))   #NO CORRELATION
cor(Jtots_nrs,deltaprimes_n2_nrs)
corspearman(Jtots_nrs,deltaprimes_n2_nrs)
savefig("/Users/sbenjamin/Desktop/Itot.png")
"""
    using PlotlyJS
    PlotlyJS.plot(histogram2dcontour(x=Jtgrad_at_rs_n2_nrs_nowellhills./1e6,y=deltaprimes_n2_nrs_nowellhills,
            contours = attr(
            showlabels = true,
            labelfont = attr(
                family = "Raleway",
                color = "white"
            )
        ),
        hoverlabel = attr(
        bgcolor = "white",
        bordercolor = "black",
        font = attr(
            family = "Raleway",
            color = "black"
        )
    )))

    trace1 = histogram2dcontour(
        x = Jtgrad_at_rs_n2_nrs_nowellhills./1e6,
        y = deltaprimes_n2_nrs_nowellhills,
        colorscale = "Blues",
        reversescale = true,
        xaxis = "x",
        yaxis = "y"
    )
    trace2 = scatter(
            x = Jtgrad_at_rs_n2_nrs_nowellhills./1e6,
            y = deltaprimes_n2_nrs_nowellhills,
            xaxis = "x",
            yaxis = "y",
            mode = "markers",
            marker = attr(
                color = "rgba(0,0,0,0.3)",
                size = 3
            )
        )
    layout = Layout(
    autosize = false,
    xaxis = attr(
        zeroline = false,
        domain = [0,1.0],
        showgrid = false
    ),
    yaxis = attr(
        zeroline = false,
        domain = [0,1.0],
        showgrid = false
    ),
    height = 600,
    width = 600,
    bargap = 0,
    hovermode = "closest",
    showlegend = false
    )

    PlotlyJS.plot([trace1, trace2], layout)
"""

plot_current_profiles(non_rs_equil_nowells[sorted_inds_nowell[1:10]];m_ind=1)
plot_current_profiles(non_rs_equil_nowellhills[sorted_inds_nowellhill[1:3]];m_ind=1)


Jts_nrs = [i.equilibrium.Jt for i in non_rs_equil]                                          #
Jps_nrs = [i.equilibrium.Jp for i in non_rs_equil]                                          #
rss_n2_nrs = [i.Δprimes[1].rs for i in non_rs_equil]                                        #
rss_n3_nrs = [i.Δprimes[2].rs for i in non_rs_equil]                                        #
Jtots_nrs = [total_plasma_current(i.equilibrium.Jt,i.equilibrium.rb) for i in non_rs_equil] #
Jt_at_rs_n2_nrs = [Jts_nrs[i](rss_n2_nrs[i]) for i in 1:length(non_rs_equil)]               
Jt_at_rs_n3_nrs = [Jts_nrs[i](rss_n3_nrs[i]) for i in 1:length(non_rs_equil)]
Jp_at_rs_n2_nrs = [Jps_nrs[i](rss_n2_nrs[i]) for i in 1:length(non_rs_equil)]
Jp_at_rs_n3_nrs = [Jps_nrs[i](rss_n3_nrs[i]) for i in 1:length(non_rs_equil)]
shear_at_rs_n2_nrs = [ForwardDiff.derivative(non_rs_equil[i].equilibrium.q,rss_n2_nrs[i]) for i in 1:length(non_rs_equil)]
shear_at_rs_n3_nrs = [ForwardDiff.derivative(non_rs_equil[i].equilibrium.q,rss_n3_nrs[i]) for i in 1:length(non_rs_equil)]
qdd_at_rs_n2_nrs = [ForwardDiff.derivative(x0->ForwardDiff.derivative(non_rs_equil[i].equilibrium.q,x0),rss_n2_nrs[i]) for i in 1:length(non_rs_equil)]
q_at_rs_n2_nrs = [non_rs_equil[i].equilibrium.q(rss_n2_nrs[i]) for i in 1:length(non_rs_equil)]
q_at_rs_n3_nrs = [non_rs_equil[i].equilibrium.q(rss_n3_nrs[i]) for i in 1:length(non_rs_equil)]
p_at_rs_n2_nrs= [mu0*non_rs_equil[i].equilibrium.p(rss_n2_nrs[i]) for i in 1:length(non_rs_equil)]
p_at_rs_n3_nrs = [mu0*non_rs_equil[i].equilibrium.p(rss_n3_nrs[i]) for i in 1:length(non_rs_equil)]
dpdr_at_rs_n2_nrs = [non_rs_equil[i].equilibrium.dpdr(rss_n2_nrs[i]) for i in 1:length(non_rs_equil)]
dpdr_at_rs_n3_nrs = [non_rs_equil[i].equilibrium.dpdr(rss_n3_nrs[i]) for i in 1:length(non_rs_equil)]



#What makes the current profile more stable... in a cylinder, its being on a Jt hill, not being in a Jt well.
    #Also the instantaneous current Jt gradient being small 
    
#Cases to confirm:

#CHOOSING EQUILIBRIA
# P1   1. A general stability to instability scan - randomly take a couple of equilibria from 5-10 segments of the delta' scale (for n=2, with n=3 stable for simplicity)
    
# P2   2. Test some interesting cases:
            #Some n=2 stable, n=3 unstable cases
            #A couple n=2, n=3 unstable cases
            #A couple n=2 marginal, n=3 wildly unstable cases
# P111 3. A sequence of well and hill equilibria (since this is my strongest argument)
# P2   4. Remake Cesar's stability to instability case 

hill_m3dc1_equils
well_m3dc1_equils

cd("/Users/sbenjamin/Desktop/PHD/ScrewPinchDataBase")
rvec = range(0.01,rb;step=0.01)
#print_equil_data(hill_m3dc1_equils; directoryprefactor="hill_equil", rvec=rvec)
#print_equil_data(well_m3dc1_equils; directoryprefactor="well_equil", rvec=rvec)

path="/Users/sbenjamin/Desktop/PHD/ScrewPinchDataBase"
for i in 1:length(hill_m3dc1_equils)
    mkdir("hill_equil_$(i)")
    cd("hill_equil_$(i)")

    print_equil_data_(hill_m3dc1_equils[i]; rvec=rvec)
    cd(path)
end

for i in 1:length(well_m3dc1_equils)
    mkdir("well_equil_$(i)")
    cd("well_equil_$(i)")

    print_equil_data_(well_m3dc1_equils[i]; rvec=rvec)
    cd(path)
end

mkdir("well_equil_1extra")
cd("well_equil_1extra")

print_equil_data_(well_m3dc1_equils_extra; rvec=rvec)
cd(path)





#make list
#make function to gen_inputs
#download the equilibria onto julia on the cluster


#RUN PLAN:
# P2   . Do some 6-field to confirm instability 






















#USE THESE FOR TITLES:
#histogram2d(grad_J_in[filtinds],deltaprimes[filtinds]; xlabel="∇J/J in",ylabel = "Δ'", title="∇J/Jin correlation with Δ'")
#cor(grad_J_in[filtinds],deltaprimes[filtinds])
#histogram2d(grad_J_out,deltaprimes; xlabel="∇J/J out",ylabel = "Δ'", title="∇J/Jout correlation with Δ'") 
#cor(grad_J_out,deltaprimes)